## In this script survey research is analyzed 

library(readxl)
install.packages("psych")
library(psych)
library(ggplot2)
## read excel and descriptive details

survey_data <- read_excel("survey_data.xlsx")
str(survey_data)
View(survey_data)

### factor analysis of numeric variables (except of social reilef)
# Specify the names of the variables you want to include in the subset
selected_vars <- c("political_int", "traditional_med", "social_med", 
                   "elite_attitude1", "elite_attitude2", "elite_attitude3", 
                   "gov_perf_earthquake1", "gov_perf_earthquake2", 
                   "pol_stab1", "pol_stab2", "pol_unstab1", "pol_stab3", 
                   "pol_stab4")

selected_vars <- setdiff(selected_vars, "social_relief")
subset_data <- survey_data[, selected_vars]
subset_data <- na.omit(subset_data)
factor_results <- fa(subset_data)
print(factor_results)

## factor analysis for social relief variables

social_data <- na.omit(social_data)
View(social_data)
# Conduct the factor analysis
factor_results <- fa(social_data, nfactors = 1, rotate = "varimax")

# Print the factor analysis results
print(factor_results)


# Check standard deviation of variables
sd_social_data <- apply(social_data, 2, sd)

# Identify variables with zero standard deviation
vars_with_zero_sd <- names(sd_social_data[sd_social_data == 0])

# Remove variables with zero standard deviation
social_data <- social_data[, !(names(social_data) %in% vars_with_zero_sd)]
View(social_data)
# Conduct the factor analysis
factor_results <- fa(social_data, nfactors = 1, rotate = "varimax")

# Print the factor analysis results
print(factor_results)

## according to results of factors, there is no strong redundancy and variability are not meaningful

View(survey_data)

## other factor analysis for other variables

selected_vars <- c("candidate_based_voting_2023", "party_based_voting_2023", "candidate_based_voting_2018", 
                   "party_based_voting_2018", "like_dislike_akp", "like_dislike_chp", 
                   "like_dislike_iyip", "like_dislike_mhp","like_dislike_hdp", "like_dislike_yrp", 
                   "like_dislike_zafp", "like_dislike_bbp", "like_dislike_other_party", "like_dislike_tkp", 
                   "like_dislike_mp", "like_dislike_rte", "like_dislike_kemal", "like_dislike_muharrem",
                   "like_dislike_sinan","satisf_democ")
selected_data <- survey_data[, selected_vars]
selected_data <- na.omit(selected_data)
factor_results <- fa(selected_data, nfactors = 1, rotate = "varimax")
print(factor_results)


### create new variables
opponents_items <- c("like_dislike_kemal", "like_dislike_muharrem", "like_dislike_sinan")
survey_data$Attitudes_to_Opponents <- rowMeans(survey_data[, opponents_items])
View(survey_data)
str(survey_data)

nation_alliance_items <- c("like_dislike_chp", "like_dislike_iyip", "like_dislike_hdp", "like_dislike_other_party")
survey_data$Attitudes_to_Nation_Alliance <- rowMeans(survey_data[, nation_alliance_items])
survey_data$Attitudes_to_Peoples_Alliance <- rowMeans(survey_data[, c("like_dislike_akp", "like_dislike_mhp", "like_dislike_yrp", "like_dislike_bbp")], na.rm = TRUE)

survey_data <- survey_data[, !colnames(survey_data) %in% c("social_relief", "social_relief_concern")]

## make again factor analysis for these variables

items <- survey_data[, c("Attitudes_to_Opponents", "Attitudes_to_Nation_Alliance")]
factor_analysis <- fa(items)
print(factor_analysis)


str(survey_data)
View(survey_data)
survey_data$age <- as.numeric(survey_data$age)

## create variables for analysis (factor analyses are made)

survey_data$government_performance <- rowMeans(survey_data[, c("gov_perf_earthquake1", "gov_perf_earthquake2")], na.rm = TRUE)
survey_data$political_stability <- rowMeans(survey_data[, c("pol_stab1", "pol_stab2")], na.rm = TRUE)


## some analysis
# Fit logistic regression model and assign it to model 1 (socio-demographic and explanatory variable)
model1 <- glm(vote_for_erd_20231 ~ age + female + education + Income + urban + Attitudes_to_Nation_Alliance+
                Attitudes_to_Peoples_Alliance+ government_performance, 
              data = survey_data, family = binomial)

# Summary of the model 1
summary(model1)

# Fit logistic regression model 2 (socio-demographic and explanatory variable)
model2 <- glm(vote_for_erd_20232 ~ age + female + education + Income + urban + Attitudes_to_Nation_Alliance+
                Attitudes_to_Peoples_Alliance+ government_performance, 
              data = survey_data, family = binomial)

# Display summary of the model2
summary(model2)
View(survey_data)
# Fit logistic regression model3 (socio-demographic, control and explanatory variable)

model3 <- glm(vote_for_erd_20231 ~ age + female + education + Income + urban + 
               government_performance + political_int + left_right_Self + earthquake_lose_person + 
                earhthquake_lose_home + Attitudes_to_Nation_Alliance + 
                Attitudes_to_Peoples_Alliance, 
              family = binomial, data = survey_data)

# Display summary of the model3
summary(model3)

# Fit logistic regression model4 (socio-demographic, conntrol and explanatory variable)

model4 <- glm(vote_for_erd_20232 ~ age + female + education + Income + urban + 
                government_performance + political_int + left_right_Self + earthquake_lose_person + 
                earhthquake_lose_home + Attitudes_to_Nation_Alliance + 
                Attitudes_to_Peoples_Alliance, 
              family = binomial, data = survey_data)


# Display summary of the model4
summary(model4)

## model 5, in this model dependent variable is vote for people's alliance
## and independent variables are mostly socio-demographic informations of respondents

## fit logistic regression model5 (socio-demographic and explanatory variable)

model5 <- glm( vote_for_peoples_alliance~ age + female + education + Income + urban + Attitudes_to_Nation_Alliance+
                 Attitudes_to_Peoples_Alliance+ government_performance, 
               data = survey_data, family = binomial)

# display the summer of model5
summary(model5)


## fit logistic regression model6 (socio-demographic, control and explanatory variable)

model6 <- glm(vote_for_peoples_alliance ~ age + female + education + Income + urban + 
                government_performance + political_int + left_right_Self + earthquake_lose_person + 
                earhthquake_lose_home + Attitudes_to_Nation_Alliance + 
                Attitudes_to_Peoples_Alliance, 
              family = binomial, data = survey_data)


# display the summer of model6
summary(model6)

### some predicted probabilities

# for model 3 (peoples alliance)
# Define the range for Attitudes_to_Peoples_Alliance
attitudes_seq <- seq(min(survey_data$Attitudes_to_Peoples_Alliance, na.rm = TRUE), 
                     max(survey_data$Attitudes_to_Peoples_Alliance, na.rm = TRUE), 
                     length.out = 100)

# Extract means for other variables
mean_values <- colMeans(survey_data[, c("age", "female", "education", "Income", "urban", 
                                        "political_int", "left_right_Self", 
                                        "earthquake_lose_person", "earhthquake_lose_home", 
                                        "Attitudes_to_Nation_Alliance")], na.rm = TRUE)

# Create a dataframe for new data
new_data <- expand.grid(Attitudes_to_Peoples_Alliance = attitudes_seq,
                        government_performance = c(min(survey_data$government_performance, na.rm = TRUE), 
                                                   max(survey_data$government_performance, na.rm = TRUE)))

# Add columns for other variables with their mean values
new_data <- cbind(new_data, 
                  age = mean_values["age"],
                  female = mean_values["female"],
                  education = mean_values["education"],
                  Income = mean_values["Income"],
                  urban = mean_values["urban"],
                  political_int = mean_values["political_int"],
                  left_right_Self = mean_values["left_right_Self"],
                  earthquake_lose_person = mean_values["earthquake_lose_person"],
                  earhthquake_lose_home = mean_values["earhthquake_lose_home"],
                  Attitudes_to_Nation_Alliance = mean_values["Attitudes_to_Nation_Alliance"])

# Predict probabilities and standard errors
preds <- predict(model3, newdata = new_data, type = "link", se.fit = TRUE)

# Compute confidence intervals
z <- qnorm(0.975)  # for 95% confidence intervals
new_data$fit <- plogis(preds$fit)  # Convert from log-odds to probability
new_data$se <- preds$se.fit
new_data$upr <- plogis(preds$fit + z * new_data$se)  # Upper bound
new_data$lwr <- plogis(preds$fit - z * new_data$se)  # Lower bound


# Create the plot with adjusted text sizes and legend position
plot <- ggplot(new_data, aes(x = Attitudes_to_Peoples_Alliance, y = fit, color = as.factor(government_performance))) +
  geom_line(size = 1) +
  geom_ribbon(aes(ymin = lwr, ymax = upr, fill = as.factor(government_performance)), alpha = 0.2, color = NA) +
  labs(title = "Attitudes to People's Alliance and Perceived Gov. Stab",
       x = "Attitudes to People's Alliance",
       y = "Pr(Vote_for_Erdoğan",
       color = "Perceived Gov. Stab.",
       fill = "Perceived Gov. Stab.") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 10, face = "bold", hjust = 0.5),  # Center the title
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 12),
    legend.title = element_text(size = 10),  # Smaller legend title
    legend.text = element_text(size = 9),    # Smaller legend text
    legend.position = "bottom",               # Move legend to the bottom
    legend.box = "horizontal",                # Arrange legend items horizontally
    legend.key.size = unit(0.5, "cm"),        # Decrease size of legend keys
    plot.margin = margin(10, 10, 30, 10)      # Adjust margins to avoid excessive blank space
  )

# Display the plot
print(plot)

# Save the plot to a file with increased dimensions
ggsave("plot_with_conf_intervals_adjusted.png", 
       plot = plot, 
       width = 12, height = 8, units = "in", 
       dpi = 300)



### predicted probability for model 3 (nation alliance)




# Generate sequence for Attitudes_to_Nation_Alliance
attitudes_seq <- seq(min(survey_data$Attitudes_to_Nation_Alliance, na.rm = TRUE), 
                     max(survey_data$Attitudes_to_Nation_Alliance, na.rm = TRUE), 
                     length.out = 100)

# Extract means for other variables
mean_values <- colMeans(survey_data[, c("age", "female", "education", "Income", "urban", 
                                        "political_int", "left_right_Self", 
                                        "earthquake_lose_person", "earhthquake_lose_home", 
                                        "Attitudes_to_Peoples_Alliance")], na.rm = TRUE)

# Create a dataframe for new data
new_data <- expand.grid(Attitudes_to_Nation_Alliance = attitudes_seq,
                        government_performance = c(min(survey_data$government_performance, na.rm = TRUE), 
                                                   max(survey_data$government_performance, na.rm = TRUE)))

# Add columns for other variables with their mean values
new_data <- cbind(new_data, 
                  age = mean_values["age"],
                  female = mean_values["female"],
                  education = mean_values["education"],
                  Income = mean_values["Income"],
                  urban = mean_values["urban"],
                  political_int = mean_values["political_int"],
                  left_right_Self = mean_values["left_right_Self"],
                  earthquake_lose_person = mean_values["earthquake_lose_person"],
                  earhthquake_lose_home = mean_values["earhthquake_lose_home"],
                  Attitudes_to_Peoples_Alliance = mean_values["Attitudes_to_Peoples_Alliance"])

# Predict probabilities and standard errors
preds <- predict(model3, newdata = new_data, type = "link", se.fit = TRUE)

# Compute confidence intervals
z <- qnorm(0.975)  # for 95% confidence intervals
new_data$fit <- plogis(preds$fit)  # Convert from log-odds to probability
new_data$se <- preds$se.fit
new_data$upr <- plogis(preds$fit + z * new_data$se)  # Upper bound
new_data$lwr <- plogis(preds$fit - z * new_data$se)  # Lower bound

# Create the plot with adjusted text sizes and legend position
plot <- ggplot(new_data, aes(x = Attitudes_to_Nation_Alliance, y = fit, color = as.factor(government_performance))) +
  geom_line(size = 1) +
  geom_ribbon(aes(ymin = lwr, ymax = upr, fill = as.factor(government_performance)), alpha = 0.2, color = NA) +
  labs(title = "Attitudes to Nation Alliance and Perceived Gov. Stab",
       x = "Attitudes to Nation Alliance",
       y = "Pr(Vote_for_Erdoğan)",
       color = "Perceived Gov. Stab.",
       fill = "Perceived Gov. Stab.") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 10, face = "bold", hjust = 0.5),  # Center the title
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 12),
    legend.title = element_text(size = 10),  # Smaller legend title
    legend.text = element_text(size = 9),    # Smaller legend text
    legend.position = "bottom",               # Move legend to the bottom
    legend.box = "horizontal",                # Arrange legend items horizontally
    legend.key.size = unit(0.5, "cm"),        # Decrease size of legend keys
    plot.margin = margin(10, 10, 30, 10)      # Adjust margins to avoid excessive blank space
  )

# Display the plot
print(plot)

# Save the plot to a file with increased dimensions
ggsave("plot_with_conf_intervals_adjusted.png", 
       plot = plot, 
       width = 12, height = 8, units = "in", 
       dpi = 300)


### predicted probability for model 6 (vote for people alliance)


library(ggplot2)

# Define the sequence for Attitudes_to_Peoples_Alliance
attitudes_seq <- seq(min(survey_data$Attitudes_to_Peoples_Alliance, na.rm = TRUE), 
                     max(survey_data$Attitudes_to_Peoples_Alliance, na.rm = TRUE), 
                     length.out = 100)

# Extract means for other variables
mean_values <- colMeans(survey_data[, c("age", "female", "education", "Income", "urban", 
                                        "political_int", "left_right_Self", 
                                        "earthquake_lose_person", "earhthquake_lose_home", 
                                        "Attitudes_to_Nation_Alliance")], na.rm = TRUE)

# Create a dataframe for new data
new_data <- expand.grid(Attitudes_to_Peoples_Alliance = attitudes_seq,
                        government_performance = c(min(survey_data$government_performance, na.rm = TRUE), 
                                                   max(survey_data$government_performance, na.rm = TRUE)))

# Add columns for other variables with their mean values
new_data <- cbind(new_data, 
                  age = mean_values["age"],
                  female = mean_values["female"],
                  education = mean_values["education"],
                  Income = mean_values["Income"],
                  urban = mean_values["urban"],
                  political_int = mean_values["political_int"],
                  left_right_Self = mean_values["left_right_Self"],
                  earthquake_lose_person = mean_values["earthquake_lose_person"],
                  earhthquake_lose_home = mean_values["earhthquake_lose_home"],
                  Attitudes_to_Nation_Alliance = mean_values["Attitudes_to_Nation_Alliance"])

# Predict probabilities and standard errors
preds <- predict(model6, newdata = new_data, type = "link", se.fit = TRUE)

# Compute confidence intervals
z <- qnorm(0.975)  # for 95% confidence intervals
new_data$fit <- plogis(preds$fit)  # Convert from log-odds to probability
new_data$se <- preds$se.fit
new_data$upr <- plogis(preds$fit + z * new_data$se)  # Upper bound
new_data$lwr <- plogis(preds$fit - z * new_data$se)  # Lower bound

# Create the plot with adjusted text sizes and legend position
plot <- ggplot(new_data, aes(x = Attitudes_to_Peoples_Alliance, y = fit, color = as.factor(government_performance))) +
  geom_line(size = 1) +
  geom_ribbon(aes(ymin = lwr, ymax = upr, fill = as.factor(government_performance)), alpha = 0.2, color = NA) +
  labs(title = "Attitudes to People's Alliance and Perceived Gov. Stab",
       x = "Attitudes to People's Alliance",
       y = "Pr(Vote_for_People's_Alliance)",
       color = "Perceived Gov. Stab.",
       fill = "Perceived Gov. Stab.") +
  
  theme_minimal() +
  theme(
    plot.title = element_text(size = 10, face = "bold", hjust = 0.5),  # Center the title
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 12),
    legend.title = element_text(size = 10),  # Smaller legend title
    legend.text = element_text(size = 9),    # Smaller legend text
    legend.position = "bottom",               # Move legend to the bottom
    legend.box = "horizontal",                # Arrange legend items horizontally
    legend.key.size = unit(0.5, "cm"),        # Decrease size of legend keys
    plot.margin = margin(10, 10, 30, 10)      # Adjust margins to avoid excessive blank space
  )

# Display the plot
print(plot)

# Save the plot to a file with increased dimensions
ggsave("plot_with_conf_intervals_adjusted_peoples_alliance.png", 
       plot = plot, 
       width = 12, height = 8, units = "in", 
       dpi = 300)


#### predicted probability for vote for people's alliance (model 6)

library(ggplot2)

# Define the sequence for Attitudes_to_Nation_Alliance
attitudes_seq <- seq(min(survey_data$Attitudes_to_Nation_Alliance, na.rm = TRUE), 
                     max(survey_data$Attitudes_to_Nation_Alliance, na.rm = TRUE), 
                     length.out = 100)

# Extract means for other variables
mean_values <- colMeans(survey_data[, c("age", "female", "education", "Income", "urban", 
                                        "political_int", "left_right_Self", 
                                        "earthquake_lose_person", "earhthquake_lose_home", 
                                        "Attitudes_to_Peoples_Alliance")], na.rm = TRUE)

# Create a dataframe for new data
new_data <- expand.grid(Attitudes_to_Nation_Alliance = attitudes_seq,
                        government_performance = c(min(survey_data$government_performance, na.rm = TRUE), 
                                                   max(survey_data$government_performance, na.rm = TRUE)))

# Add columns for other variables with their mean values
new_data <- cbind(new_data, 
                  age = mean_values["age"],
                  female = mean_values["female"],
                  education = mean_values["education"],
                  Income = mean_values["Income"],
                  urban = mean_values["urban"],
                  political_int = mean_values["political_int"],
                  left_right_Self = mean_values["left_right_Self"],
                  earthquake_lose_person = mean_values["earthquake_lose_person"],
                  earhthquake_lose_home = mean_values["earhthquake_lose_home"],
                  Attitudes_to_Peoples_Alliance = mean_values["Attitudes_to_Peoples_Alliance"])

# Predict probabilities and standard errors
preds <- predict(model6, newdata = new_data, type = "link", se.fit = TRUE)

# Compute confidence intervals
z <- qnorm(0.975)  # for 95% confidence intervals
new_data$fit <- plogis(preds$fit)  # Convert from log-odds to probability
new_data$se <- preds$se.fit
new_data$upr <- plogis(preds$fit + z * new_data$se)  # Upper bound
new_data$lwr <- plogis(preds$fit - z * new_data$se)  # Lower bound

# Create the plot with adjusted text sizes and legend position
plot <- ggplot(new_data, aes(x = Attitudes_to_Nation_Alliance, y = fit, color = as.factor(government_performance))) +
  geom_line(size = 1) +
  geom_ribbon(aes(ymin = lwr, ymax = upr, fill = as.factor(government_performance)), alpha = 0.2, color = NA) +
  labs(title = "Attitudes to Nation Alliance and Perceived Gov. Stab.",
       x = "Attitudes to Nation Alliance",
       y = "Pr(Vote_for_People's_Alliance)",
       color = "Perceived Gov. Stab.",
       fill = "Perceived Gov. Stab.") +
  
  theme_minimal() +
  theme(
    plot.title = element_text(size = 10, face = "bold", hjust = 0.5),  # Center the title
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 12),
    legend.title = element_text(size = 10),  # Smaller legend title
    legend.text = element_text(size = 9),    # Smaller legend text
    legend.position = "bottom",               # Move legend to the bottom
    legend.box = "horizontal",                # Arrange legend items horizontally
    legend.key.size = unit(0.5, "cm"),        # Decrease size of legend keys
    plot.margin = margin(10, 10, 30, 10)      # Adjust margins to avoid excessive blank space
  )

# Display the plot
print(plot)

# Save the plot to a file with increased dimensions
ggsave("plot_with_conf_intervals_adjusted_nation_alliance.png", 
       plot = plot, 
       width = 12, height = 8, units = "in", 
       dpi = 300)



















